      subroutine maxwell(nvtx, nvtxkm, ijkvtx, iwid, jwid, kwid, precon, tiny, &
         ncells, ijkcell, itdim, nsp, ibp1, jbp1, kbp1, vvolaxis, divetil, cdlt&
         , sdlt, strait, dz, clite, cntr, dt, fourpi, itmax, error, periodicx, &
         c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, &
         c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vvol, vol, qom, qdnv&
         , curlx, curly, curlz, exmu, eymu, ezmu, rx, ry, rz, qx, qy, qz, qxm1&
         , qym1, qzm1, aqxm1, aqym1, aqzm1, aqx, aqy, aqz, paqx, paqy, paqz, &
         a11, a12, a13, a21, a22, a23, a31, a32, a33, ijkvtmp, qxtilde, qytilde&
         , qztilde, aqxtilde, aqytilde, aqztilde, jxtilde, jytilde, jztilde, &
         paqxm1, paqym1, paqzm1, gdivx, gdivy, gdivz, qxm2, qym2, qzm2, aqxm2, &
         aqym2, aqzm2, paqxm2, paqym2, paqzm2, qxm3, qym3, qzm3, aqxm3, aqym3, &
         aqzm3, paqxm3, paqym3, paqzm3, bxv, byv, bzv, ex0, ey0, ez0, ex, ey, &
         ez, iter) 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: nvtx 
      integer  :: nvtxkm 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      integer  :: ncells 
      integer  :: itdim 
      integer  :: nsp 
      integer  :: ibp1 
      integer  :: jbp1 
      integer  :: kbp1 
      integer , intent(in) :: itmax 
      integer , intent(out) :: iter 
      real  :: tiny 
      real  :: cdlt 
      real  :: sdlt 
      real  :: strait 
      real  :: dz 
      real  :: clite 
      real  :: cntr 
      real  :: dt 
      real  :: fourpi 
      real , intent(in) :: error 
      logical , intent(in) :: precon 
      logical  :: periodicx 
      integer  :: ijkvtx(*) 
      integer  :: ijkcell(*) 
      integer  :: ijkvtmp(*) 
      real  :: vvolaxis(*) 
      real  :: divetil(*) 
      real  :: c1x(*) 
      real  :: c2x(*) 
      real  :: c3x(*) 
      real  :: c4x(*) 
      real  :: c5x(*) 
      real  :: c6x(*) 
      real  :: c7x(*) 
      real  :: c8x(*) 
      real  :: c1y(*) 
      real  :: c2y(*) 
      real  :: c3y(*) 
      real  :: c4y(*) 
      real  :: c5y(*) 
      real  :: c6y(*) 
      real  :: c7y(*) 
      real  :: c8y(*) 
      real  :: c1z(*) 
      real  :: c2z(*) 
      real  :: c3z(*) 
      real  :: c4z(*) 
      real  :: c5z(*) 
      real  :: c6z(*) 
      real  :: c7z(*) 
      real  :: c8z(*) 
      real  :: vvol(*) 
      real  :: vol(*) 
      real  :: qom(*) 
      real  :: qdnv(itdim,*) 
      real  :: curlx(*) 
      real  :: curly(*) 
      real  :: curlz(*) 
      real  :: exmu(*) 
      real  :: eymu(*) 
      real  :: ezmu(*) 
      real  :: rx(*) 
      real  :: ry(*) 
      real  :: rz(*) 
      real , intent(inout) :: qx(*) 
      real , intent(inout) :: qy(*) 
      real , intent(inout) :: qz(*) 
      real , intent(inout) :: qxm1(*) 
      real , intent(inout) :: qym1(*) 
      real , intent(inout) :: qzm1(*) 
      real , intent(inout) :: aqxm1(*) 
      real , intent(inout) :: aqym1(*) 
      real , intent(inout) :: aqzm1(*) 
      real  :: aqx(*) 
      real  :: aqy(*) 
      real  :: aqz(*) 
      real  :: paqx(*) 
      real  :: paqy(*) 
      real  :: paqz(*) 
      real  :: a11(*) 
      real  :: a12(*) 
      real  :: a13(*) 
      real  :: a21(*) 
      real  :: a22(*) 
      real  :: a23(*) 
      real  :: a31(*) 
      real  :: a32(*) 
      real  :: a33(*) 
      real  :: qxtilde(*) 
      real  :: qytilde(*) 
      real  :: qztilde(*) 
      real  :: aqxtilde(*) 
      real  :: aqytilde(*) 
      real  :: aqztilde(*) 
      real  :: jxtilde(*) 
      real  :: jytilde(*) 
      real  :: jztilde(*) 
      real , intent(inout) :: paqxm1(*) 
      real , intent(inout) :: paqym1(*) 
      real , intent(inout) :: paqzm1(*) 
      real , intent(in) :: gdivx(*) 
      real , intent(in) :: gdivy(*) 
      real , intent(in) :: gdivz(*) 
      real , intent(inout) :: qxm2(*) 
      real , intent(inout) :: qym2(*) 
      real , intent(inout) :: qzm2(*) 
      real , intent(inout) :: aqxm2(*) 
      real , intent(inout) :: aqym2(*) 
      real , intent(inout) :: aqzm2(*) 
      real , intent(inout) :: paqxm2(*) 
      real , intent(inout) :: paqym2(*) 
      real , intent(inout) :: paqzm2(*) 
      real , intent(inout) :: qxm3(*) 
      real , intent(inout) :: qym3(*) 
      real , intent(inout) :: qzm3(*) 
      real , intent(inout) :: aqxm3(*) 
      real , intent(inout) :: aqym3(*) 
      real , intent(inout) :: aqzm3(*) 
      real , intent(inout) :: paqxm3(*) 
      real , intent(inout) :: paqym3(*) 
      real , intent(inout) :: paqzm3(*) 
      real  :: bxv(*) 
      real  :: byv(*) 
      real  :: bzv(*) 
      real  :: ex0(*) 
      real  :: ey0(*) 
      real  :: ez0(*) 
      real  :: ex(*) 
      real  :: ey(*) 
      real  :: ez(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: ibp2, jbp2, kbp2, n, ijk, k, j, i, k2, ijk2, i2, nvtmp, nstart 
      real :: jx0, jy0, jz0, lambda, lambm1, jxtildet, jytildet, jztildet, &
         lambm2, lambm3, dnorm, rsum, bottom, botm1, botm2, botm3, dero, dern, &
         derd, wsqx, wsqy, wsqz, wsaqx, wsaqy, wsaqz, rnorm, alfa, dxnorm, &
         xnorm 
      logical :: testgl, precongl 
!-----------------------------------------------
!
!     solves equation for E field derived by combining Ampere's
!     law and Faraday's law with the momentum equations
!
!     csq*(cntr*dt)**2*curl(curl(E(dt+cntr*dt))) + (1.+mu)E(t+cntr*dt)
!                                = E(t) + c*cntr*dt*(curl(B(t))-jtilde)
!     MAXWELL calls AVEC
!
!
!
!
!
      testgl = .F. 
      precongl = .F. 
      if (precongl) testgl = .T. 
!
      ibp2 = ibp1 + 1 
      jbp2 = jbp1 + 1 
      kbp2 = kbp1 + 1 
!
      call axisvol (ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, ijkvtx, vol, &
         vvolaxis) 
!
!     calculate A*E0
!
!     ****************************************************************
!
!     CONSTRUCT JACOBI PRECONDITIONER
!
!     ***************************************************************
!
!
      if (precon) then 
!
         call block (nvtx, ijkvtx, nsp, itdim, iwid, jwid, kwid, c1x, c2x, c3x&
            , c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, &
            c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vvolaxis, vol, vvol, clite&
            , dt, qom, qdnv, cntr, bxv, byv, bzv, exmu, eymu, ezmu, a11, a12, &
            a13, a21, a22, a23, a31, a32, a33) 
!
!
         call ludecomp (nvtx, ijkvtx, a11, a12, a13, a21, a22, a23, a31, a32, &
            a33) 
!
      else 
!
!     write(*,*)'normal'
!        read(*,*)dnorm
         dnorm = 1. 
!dir$ ivdep
         do n = 1, nvtx 
!
            ijk = ijkvtx(n) 
!
            a11(ijk) = dnorm 
            a12(ijk) = 0. 
            a13(ijk) = 0. 
!cc
            a21(ijk) = 0. 
            a22(ijk) = dnorm 
            a23(ijk) = 0. 
!ccc
            a31(ijk) = 0. 
            a32(ijk) = 0. 
            a33(ijk) = dnorm 
!
         end do 
!
      endif 
!
!     END PRE-CONDITIIONER
!
!     ********************************************************************
!
!     CALCULATE INITIAL RESIDUE
!
!     **********************************************************************
!
      call avec (nvtx, nvtxkm, ijkvtx, iwid, jwid, kwid, tiny, ncells, ijkcell&
         , qom, qdnv, nsp, itdim, vvol, ibp1, jbp1, kbp1, vvolaxis, cdlt, sdlt&
         , strait, dz, periodicx, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, &
         c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
         c8z, vol, curlx, curly, curlz, exmu, eymu, ezmu, bxv, byv, bzv, cntr, &
         dt, 1., clite, fourpi, ex, ey, ez, rx, ry, rz) 
!
!dir$ ivdep
      do n = 1, nvtx 
!
         ijk = ijkvtx(n) 
!
!
         rx(ijk) = gdivx(ijk) - rx(ijk) 
         ry(ijk) = gdivy(ijk) - ry(ijk) 
         rz(ijk) = gdivz(ijk) - rz(ijk) 
!
      end do 
!
      call axisvec (ibp1, jbp1, iwid, jwid, kwid, vvol, vvolaxis, rx, ry, rz) 
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         rx, ry, rz, periodicx) 
!
!
!     set residue to zero on k=kbp2 plane
!
      k = kbp1 + 1 
      do j = 1, jbp1 + 1 
!
         rx(1+(j-1)*jwid+(k-1)*kwid:ibp1*iwid+1+(j-1)*jwid+(k-1)*kwid:iwid) = &
            0.0 
         ry(1+(j-1)*jwid+(k-1)*kwid:ibp1*iwid+1+(j-1)*jwid+(k-1)*kwid:iwid) = &
            0.0 
         rz(1+(j-1)*jwid+(k-1)*kwid:ibp1*iwid+1+(j-1)*jwid+(k-1)*kwid:iwid) = &
            0.0 
      end do 
!
      k = 1 
      k2 = 2 
      do j = 1, jbp1 + 1 
!
         rx(1+(j-1)*jwid+(k-1)*kwid:ibp1*iwid+1+(j-1)*jwid+(k-1)*kwid:iwid) = &
            0.0 
         ry(1+(j-1)*jwid+(k-1)*kwid:ibp1*iwid+1+(j-1)*jwid+(k-1)*kwid:iwid) = &
            0.0 
         rz(1+(j-1)*jwid+(k-1)*kwid:ibp1*iwid+1+(j-1)*jwid+(k-1)*kwid:iwid) = &
            0.0 
         rx(1+(j-1)*jwid+(k2-1)*kwid:ibp1*iwid+1+(j-1)*jwid+(k2-1)*kwid:iwid)&
             = 0.0 
         ry(1+(j-1)*jwid+(k2-1)*kwid:ibp1*iwid+1+(j-1)*jwid+(k2-1)*kwid:iwid)&
             = 0.0 
         rz(1+(j-1)*jwid+(k2-1)*kwid:ibp1*iwid+1+(j-1)*jwid+(k2-1)*kwid:iwid)&
             = 0.0 
      end do 
!
!
!     set residue to zero on i=ibp2 and i=1 plane
!
      if (.not.periodicx) then 
!
         i = ibp1 + 1 
         do j = 1, jbp1 + 1 
!
            rx(1+(i-1)*iwid+(j-1)*jwid:kbp1*kwid+1+(i-1)*iwid+(j-1)*jwid:kwid)&
                = 0.0 
            ry(1+(i-1)*iwid+(j-1)*jwid:kbp1*kwid+1+(i-1)*iwid+(j-1)*jwid:kwid)&
                = 0.0 
            rz(1+(i-1)*iwid+(j-1)*jwid:kbp1*kwid+1+(i-1)*iwid+(j-1)*jwid:kwid)&
                = 0.0 
!      ey(ijk)=0.0
!      ez(ijk)=0.0
         end do 
!
         i = 1 
         i2 = 2 
         do j = 1, jbp1 + 1 
!
            rx(1+(i-1)*iwid+(j-1)*jwid:kbp1*kwid+1+(i-1)*iwid+(j-1)*jwid:kwid)&
                = 0.0 
            ry(1+(i-1)*iwid+(j-1)*jwid:kbp1*kwid+1+(i-1)*iwid+(j-1)*jwid:kwid)&
                = 0.0 
            rz(1+(i-1)*iwid+(j-1)*jwid:kbp1*kwid+1+(i-1)*iwid+(j-1)*jwid:kwid)&
                = 0.0 
!      ey(ijk)=0.0
!      ez(ijk)=0.0
            rx(1+(i2-1)*iwid+(j-1)*jwid:kbp1*kwid+1+(i2-1)*iwid+(j-1)*jwid:kwid&
               ) = 0.0 
            ry(1+(i2-1)*iwid+(j-1)*jwid:kbp1*kwid+1+(i2-1)*iwid+(j-1)*jwid:kwid&
               ) = 0.0 
            rz(1+(i2-1)*iwid+(j-1)*jwid:kbp1*kwid+1+(i2-1)*iwid+(j-1)*jwid:kwid&
               ) = 0.0 
!      ey(ijk2)=0.0
!      ez(ijk2)=0.0
         end do 
!
      endif 
 
!
!
!     END CALCULATION OF RESIDUE
!
!     ****************************************************************
!
      if (periodicx) then 
         call list (2, ibp1, 2, jbp1, 3, kbp1, iwid, jwid, kwid, nvtmp, ijkvtmp&
            ) 
      else 
         call list (3, ibp1, 2, jbp1, 3, kbp1, iwid, jwid, kwid, nvtmp, ijkvtmp&
            ) 
      endif 
!
      qxtilde(:ibp2*jbp2*kbp2) = ex(:ibp2*jbp2*kbp2) 
      qytilde(:ibp2*jbp2*kbp2) = ey(:ibp2*jbp2*kbp2) 
      qztilde(:ibp2*jbp2*kbp2) = ez(:ibp2*jbp2*kbp2) 
!
!
      rsum = 0.0 
! ijkgmx=0
!dir$ ivdep
      do n = 1, nvtmp 
!
         ijk = ijkvtmp(n) 
!
         qx(ijk) = 0.0 
         aqx(ijk) = 0.0 
         paqx(ijk) = 0.0 
!
         qy(ijk) = 0.0 
         aqy(ijk) = 0.0 
         paqy(ijk) = 0.0 
!
         qz(ijk) = 0.0 
         aqz(ijk) = 0.0 
         paqz(ijk) = 0.0 
!
         qxm1(ijk) = 0.0 
         aqxm1(ijk) = 0.0 
         paqxm1(ijk) = 0.0 
!
         qym1(ijk) = 0.0 
         aqym1(ijk) = 0.0 
         paqym1(ijk) = 0.0 
!
         qzm1(ijk) = 0.0 
         aqzm1(ijk) = 0.0 
         paqzm1(ijk) = 0.0 
!
         qxm2(ijk) = 0.0 
         aqxm2(ijk) = 0.0 
         paqxm2(ijk) = 0.0 
!
         qym2(ijk) = 0.0 
         aqym2(ijk) = 0.0 
         paqym2(ijk) = 0.0 
!
         qzm2(ijk) = 0.0 
         aqzm2(ijk) = 0.0 
         paqzm2(ijk) = 0.0 
!
         qxm3(ijk) = 0.0 
         aqxm3(ijk) = 0.0 
         paqxm3(ijk) = 0.0 
!
         qym3(ijk) = 0.0 
         aqym3(ijk) = 0.0 
         paqym3(ijk) = 0.0 
!
         qzm3(ijk) = 0.0 
         aqzm3(ijk) = 0.0 
         paqzm3(ijk) = 0.0 
!
!
         rsum = rsum + abs(gdivx(ijk)) + abs(gdivy(ijk)) + abs(gdivz(ijk)) 
!
!
!
      end do 
!
      if (rsum == 0.0) then 
         iter = 0 
         go to 10 
      endif 
!
      bottom = 0.0 
      botm1 = 0.0 
      botm2 = 0.0 
      botm3 = 0.0 
!
!
!     **************  MAIN ITERATION LOOP  **************************
!     begin conjugate gradient interation
!
      nstart = 1 
!
      do iter = 1, itmax 
!
!
!     APPLY BLOCK PRE-CONDITIONING
!
!
         call precond (nvtmp, ijkvtmp, a11, a12, a13, a21, a22, a23, a31, a32, &
            a33, rx, ry, rz, qxtilde, qytilde, qztilde) 
!
         call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, &
            cdlt, qxtilde, qytilde, qztilde, periodicx) 
!
!
!
!dir$ ivdep
         do n = nstart, nvtmp 
!
            ijk = ijkvtmp(n) 
!
            qxtilde(ijk) = ex(ijk) + qxtilde(ijk) 
            qytilde(ijk) = ey(ijk) + qytilde(ijk) 
            qztilde(ijk) = ez(ijk) + qztilde(ijk) 
!
         end do 
!
         call axisvec (ibp1, jbp1, iwid, jwid, kwid, vvol, vvolaxis, qxtilde, &
            qytilde, qztilde) 
         call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, &
            cdlt, qxtilde, qytilde, qztilde, periodicx) 
!
         call avec (nvtx, nvtxkm, ijkvtx, iwid, jwid, kwid, tiny, nvtmp, &
            ijkvtmp, qom, qdnv, nsp, itdim, vvol, ibp1, jbp1, kbp1, vvolaxis, &
            cdlt, sdlt, strait, dz, periodicx, c1x, c2x, c3x, c4x, c5x, c6x, &
            c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
            c4z, c5z, c6z, c7z, c8z, vol, curlx, curly, curlz, exmu, eymu, ezmu&
            , bxv, byv, bzv, cntr, dt, 1., clite, fourpi, qxtilde, qytilde, &
            qztilde, aqxtilde, aqytilde, aqztilde) 
!
!dir$ ivdep
         if (precongl) then 
            do n = nstart, nvtmp 
!
               ijk = ijkvtmp(n) 
!
               aqxtilde(ijk) = rx(ijk) - (gdivx(ijk)-aqxtilde(ijk)) 
               aqytilde(ijk) = ry(ijk) - (gdivy(ijk)-aqytilde(ijk)) 
               aqztilde(ijk) = rz(ijk) - (gdivz(ijk)-aqztilde(ijk)) 
!
               ex(ijk) = qxtilde(ijk) 
               ey(ijk) = qytilde(ijk) 
               ez(ijk) = qztilde(ijk) 
!
               qxtilde(ijk) = qxtilde(ijk) - ex(ijk) 
               qytilde(ijk) = qytilde(ijk) - ey(ijk) 
               qztilde(ijk) = qztilde(ijk) - ez(ijk) 
!
            end do 
         else 
            do n = nstart, nvtmp 
!
               ijk = ijkvtmp(n) 
!
               aqxtilde(ijk) = rx(ijk) - (gdivx(ijk)-aqxtilde(ijk)) 
               aqytilde(ijk) = ry(ijk) - (gdivy(ijk)-aqytilde(ijk)) 
               aqztilde(ijk) = rz(ijk) - (gdivz(ijk)-aqztilde(ijk)) 
!
               qxtilde(ijk) = qxtilde(ijk) - ex(ijk) 
               qytilde(ijk) = qytilde(ijk) - ey(ijk) 
               qztilde(ijk) = qztilde(ijk) - ez(ijk) 
!
            end do 
         endif 
!
!     Test put on by GL & JUB
!
         dero = 0. 
         dern = 0. 
         derd = 0. 
         do n = nstart, nvtmp 
!
            ijk = ijkvtmp(n) 
!
            derd = derd + abs(aqxtilde(ijk)) + abs(aqytilde(ijk)) + abs(&
               aqztilde(ijk)) 
!
            dero = dero + abs(rx(ijk)) + abs(ry(ijk)) + abs(rz(ijk)) 
!
            dern = dern + abs(aqxtilde(ijk)-rx(ijk)) + abs(aqytilde(ijk)-ry(ijk&
               )) + abs(aqztilde(ijk)-rz(ijk)) 
!
         end do 
!      write(*,*)'TEST Conv:',iter,'r=',dero,'rn=',dern
!      write(*,*)'TEST Conv:','aq=',derd
!
         lambda = 0.0 
         lambm1 = 0.0 
         lambm2 = 0.0 
         lambm3 = 0.0 
!
!dir$ ivdep
         do n = nstart, nvtmp 
!
            ijk = ijkvtmp(n) 
!
            lambda = lambda + aqxtilde(ijk)*paqx(ijk) + aqytilde(ijk)*paqy(ijk)&
                + aqztilde(ijk)*paqz(ijk) 
!
            lambm1 = lambm1 + aqxtilde(ijk)*paqxm1(ijk) + aqytilde(ijk)*paqym1(&
               ijk) + aqztilde(ijk)*paqzm1(ijk) 
!
            lambm2 = lambm2 + aqxtilde(ijk)*paqxm2(ijk) + aqytilde(ijk)*paqym2(&
               ijk) + aqztilde(ijk)*paqzm2(ijk) 
!
            lambm3 = lambm3 + aqxtilde(ijk)*paqxm3(ijk) + aqytilde(ijk)*paqym3(&
               ijk) + aqztilde(ijk)*paqzm3(ijk) 
!
!
         end do 
!
         lambda = lambda/(bottom + 1.E-100) 
         lambm1 = lambm1/(botm1 + 1.E-100) 
         lambm2 = lambm2/(botm2 + 1.E-200) 
         lambm3 = lambm3/(botm3 + 1.E-300) 
!
!dir$ ivdep
         do n = nstart, nvtmp 
!
            ijk = ijkvtmp(n) 
!
!
            wsqx = qxtilde(ijk) - lambda*qx(ijk) - lambm1*qxm1(ijk) - lambm2*&
               qxm2(ijk) - lambm3*qxm3(ijk) 
            wsqy = qytilde(ijk) - lambda*qy(ijk) - lambm1*qym1(ijk) - lambm2*&
               qym2(ijk) - lambm3*qym3(ijk) 
            wsqz = qztilde(ijk) - lambda*qz(ijk) - lambm1*qzm1(ijk) - lambm2*&
               qzm2(ijk) - lambm3*qzm3(ijk) 
!
            wsaqx = aqxtilde(ijk) - lambda*aqx(ijk) - lambm1*aqxm1(ijk) - &
               lambm2*aqxm2(ijk) - lambm3*aqxm3(ijk) 
            wsaqy = aqytilde(ijk) - lambda*aqy(ijk) - lambm1*aqym1(ijk) - &
               lambm2*aqym2(ijk) - lambm3*aqym3(ijk) 
            wsaqz = aqztilde(ijk) - lambda*aqz(ijk) - lambm1*aqzm1(ijk) - &
               lambm2*aqzm2(ijk) - lambm3*aqzm3(ijk) 
!
            qxm3(ijk) = qxm2(ijk) 
            aqxm3(ijk) = aqxm2(ijk) 
            paqxm3(ijk) = paqxm2(ijk) 
!
            qym3(ijk) = qym2(ijk) 
            aqym3(ijk) = aqym2(ijk) 
            paqym3(ijk) = paqym2(ijk) 
!
            qzm3(ijk) = qzm2(ijk) 
            aqzm3(ijk) = aqzm2(ijk) 
            paqzm3(ijk) = paqzm2(ijk) 
!
            qxm2(ijk) = qxm1(ijk) 
            aqxm2(ijk) = aqxm1(ijk) 
            paqxm2(ijk) = paqxm1(ijk) 
!
            qym2(ijk) = qym1(ijk) 
            aqym2(ijk) = aqym1(ijk) 
            paqym2(ijk) = paqym1(ijk) 
!
            qzm2(ijk) = qzm1(ijk) 
            aqzm2(ijk) = aqzm1(ijk) 
            paqzm2(ijk) = paqzm1(ijk) 
!
            qxm1(ijk) = qx(ijk) 
            aqxm1(ijk) = aqx(ijk) 
            paqxm1(ijk) = paqx(ijk) 
!
            qym1(ijk) = qy(ijk) 
            aqym1(ijk) = aqy(ijk) 
            paqym1(ijk) = paqy(ijk) 
!
            qzm1(ijk) = qz(ijk) 
            aqzm1(ijk) = aqz(ijk) 
            paqzm1(ijk) = paqz(ijk) 
!
            qx(ijk) = wsqx 
            aqx(ijk) = wsaqx 
!
            qy(ijk) = wsqy 
            aqy(ijk) = wsaqy 
!
            qz(ijk) = wsqz 
            aqz(ijk) = wsaqz 
!
         end do 
!
         call precond (nvtmp, ijkvtmp, a11, a12, a13, a21, a22, a23, a31, a32, &
            a33, aqx, aqy, aqz, paqx, paqy, paqz) 
!
         call axisvec (ibp1, jbp1, iwid, jwid, kwid, vvol, vvolaxis, paqx, paqy&
            , paqz) 
!
!
!
         rnorm = 0.0 
         alfa = 0.0 
         bottom = 0.0 
         botm1 = 0.0 
         botm2 = 0.0 
         botm3 = 0.0 
!
!
!
!dir$ ivdep
         do n = nstart, nvtmp 
!
            ijk = ijkvtmp(n) 
!
!
            rnorm = rnorm + abs(rx(ijk)) + abs(ry(ijk)) + abs(rz(ijk)) 
!
            alfa = alfa + rx(ijk)*paqx(ijk) + ry(ijk)*paqy(ijk) + rz(ijk)*paqz(&
               ijk) 
!
            bottom = bottom + aqx(ijk)*paqx(ijk) + aqy(ijk)*paqy(ijk) + aqz(ijk&
               )*paqz(ijk) 
!
            botm1 = botm1 + aqxm1(ijk)*paqxm1(ijk) + aqym1(ijk)*paqym1(ijk) + &
               aqzm1(ijk)*paqzm1(ijk) 
!
            botm2 = botm2 + aqxm2(ijk)*paqxm2(ijk) + aqym2(ijk)*paqym2(ijk) + &
               aqzm2(ijk)*paqzm2(ijk) 
!
            botm3 = botm3 + aqxm3(ijk)*paqxm3(ijk) + aqym3(ijk)*paqym3(ijk) + &
               aqzm3(ijk)*paqzm3(ijk) 
!
!
         end do 
!
         alfa = alfa/(bottom + 1.E-100) 
!
!dir$ ivdep
         if (.not.precongl) then 
            do n = nstart, nvtmp 
!
               ijk = ijkvtmp(n) 
!
               ex(ijk) = ex(ijk) + alfa*qx(ijk) 
               ey(ijk) = ey(ijk) + alfa*qy(ijk) 
               ez(ijk) = ez(ijk) + alfa*qz(ijk) 
!
               rx(ijk) = rx(ijk) - alfa*aqx(ijk) 
               ry(ijk) = ry(ijk) - alfa*aqy(ijk) 
               rz(ijk) = rz(ijk) - alfa*aqz(ijk) 
!
            end do 
         else 
            do n = nstart, nvtmp 
!
               ijk = ijkvtmp(n) 
!
               rx(ijk) = rx(ijk) - alfa*aqx(ijk) 
               ry(ijk) = ry(ijk) - alfa*aqy(ijk) 
               rz(ijk) = rz(ijk) - alfa*aqz(ijk) 
!
            end do 
         endif 
!
!     Added by GL
!
         if (testgl) then 
!
            rnorm = 0.0 
!
!dir$ ivdep
            do n = nstart, nvtmp 
!
               ijk = ijkvtmp(n) 
!
               rnorm = rnorm + abs(rx(ijk)) + abs(ry(ijk)) + abs(rz(ijk)) 
!
            end do 
            write (*, *) 'Conv, pre, rnorm=', rnorm 
!
!
            call axisvec (ibp1, jbp1, iwid, jwid, kwid, vvol, vvolaxis, ex, ey&
               , ez) 
!
            call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, &
               cdlt, ex, ey, ez, periodicx) 
!
            call avec (nvtx, nvtxkm, ijkvtx, iwid, jwid, kwid, tiny, ncells, &
               ijkcell, qom, qdnv, nsp, itdim, vvol, ibp1, jbp1, kbp1, vvolaxis&
               , cdlt, sdlt, strait, dz, periodicx, c1x, c2x, c3x, c4x, c5x, &
               c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z&
               , c3z, c4z, c5z, c6z, c7z, c8z, vol, curlx, curly, curlz, exmu, &
               eymu, ezmu, bxv, byv, bzv, cntr, dt, 1., clite, fourpi, ex, ey, &
               ez, rx, ry, rz) 
!
!dir$ ivdep
            do n = 1, nvtx 
!
               ijk = ijkvtx(n) 
!
!
               rx(ijk) = gdivx(ijk) - rx(ijk) 
               ry(ijk) = gdivy(ijk) - ry(ijk) 
               rz(ijk) = gdivz(ijk) - rz(ijk) 
!
            end do 
         endif 
!
!     END added by GL
!
!
         dxnorm = 0.0 
         xnorm = 0.0 
         if (testgl) rnorm = 0.0 
!
!dir$ ivdep
         if (testgl) then 
            do n = nstart, nvtmp 
!
               ijk = ijkvtmp(n) 
!
               rnorm = rnorm + abs(rx(ijk)) + abs(ry(ijk)) + abs(rz(ijk)) 
!
               dxnorm = dxnorm + abs(qx(ijk)) + abs(qy(ijk)) + abs(qz(ijk)) 
!
               xnorm = xnorm + abs(ex(ijk)) + abs(ey(ijk)) + abs(ez(ijk)) 
!
            end do 
         else 
            do n = nstart, nvtmp 
!
               ijk = ijkvtmp(n) 
!
               dxnorm = dxnorm + abs(qx(ijk)) + abs(qy(ijk)) + abs(qz(ijk)) 
!
               xnorm = xnorm + abs(ex(ijk)) + abs(ey(ijk)) + abs(ez(ijk)) 
!
            end do 
         endif 
!
         dxnorm = dxnorm*abs(alfa) 
!
!      write(*,*)'TEST Conv2','dxnorm=',dxnorm,'xnorm=',xnorm
!      write(*,*)'TEST Conv2','rnorm=',rnorm,'rsum=',rsum,error
         if (dxnorm<=error*xnorm .and. rnorm<=error*rsum) go to 10 
!
         call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, &
            cdlt, rx, ry, rz, periodicx) 
!
!
!
      end do 
!
      write (*, *) 'maxwell did not converge', rnorm 
!      ex(ijkvtx(:nvtx)) = 0. 
!      ey(ijkvtx(:nvtx)) = 0. 
!      ez(ijkvtx(:nvtx)) = 0. 
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         ex, ey, ez, periodicx) 
!
!
      return  
!     **********************End of Iteration **************************
!
   10 continue 
      write (*, *) 'maxwell converged', iter, rnorm 
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         ex, ey, ez, periodicx) 
!
      call axisvec (ibp1, jbp1, iwid, jwid, kwid, vvol, vvolaxis, ex, ey, ez) 
!
!
!     ********************** END MAXWELL **********************
!
      return  
      end subroutine maxwell 


      subroutine avec(nvtx, nvtxkm, ijkvtx, iwid, jwid, kwid, tiny, ncells, &
         ijkcell, qom, qdnv, nsp, itdim, vvol, ibp1, jbp1, kbp1, vvolaxis, cdlt&
         , sdlt, strait, dz, periodicx, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x&
         , c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z&
         , c7z, c8z, vol, curlvx, curlvy, curlvz, exmu, eymu, ezmu, bxv, byv, &
         bzv, cntr, dt, transpos, clite, fourpi, ex, ey, ez, hx, hy, hz) 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: nvtx 
      integer  :: nvtxkm 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      integer  :: ncells 
      integer  :: nsp 
      integer  :: itdim 
      integer  :: ibp1 
      integer  :: jbp1 
      integer  :: kbp1 
      real  :: tiny 
      real  :: cdlt 
      real  :: sdlt 
      real  :: strait 
      real  :: dz 
      real , intent(in) :: cntr 
      real  :: dt 
      real  :: transpos 
      real  :: clite 
      real  :: fourpi 
      logical  :: periodicx 
      integer  :: ijkvtx(*) 
      integer  :: ijkcell(*) 
      real  :: qom(*) 
      real  :: qdnv(itdim,*) 
      real , intent(in) :: vvol(*) 
      real , intent(in) :: vvolaxis(*) 
      real  :: c1x(*) 
      real  :: c2x(*) 
      real  :: c3x(*) 
      real  :: c4x(*) 
      real  :: c5x(*) 
      real  :: c6x(*) 
      real  :: c7x(*) 
      real  :: c8x(*) 
      real  :: c1y(*) 
      real  :: c2y(*) 
      real  :: c3y(*) 
      real  :: c4y(*) 
      real  :: c5y(*) 
      real  :: c6y(*) 
      real  :: c7y(*) 
      real  :: c8y(*) 
      real  :: c1z(*) 
      real  :: c2z(*) 
      real  :: c3z(*) 
      real  :: c4z(*) 
      real  :: c5z(*) 
      real  :: c6z(*) 
      real  :: c7z(*) 
      real  :: c8z(*) 
      real  :: vol(*) 
      real  :: curlvx(*) 
      real  :: curlvy(*) 
      real  :: curlvz(*) 
      real  :: exmu(*) 
      real  :: eymu(*) 
      real  :: ezmu(*) 
      real  :: bxv(*) 
      real  :: byv(*) 
      real  :: bzv(*) 
      real  :: ex(*) 
      real  :: ey(*) 
      real  :: ez(*) 
      real , intent(out) :: hx(*) 
      real , intent(out) :: hy(*) 
      real , intent(out) :: hz(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: j, i, ijk, n 
      real :: rvvol, dtsq, dtc, dtcsq 
      logical :: cartesian 
!-----------------------------------------------
!***BEGIN PROLOGUE
!***PURPOSE
!
!     calculate A(transpose)*A*E
!
!s***DESCRIPTION
!
!     on entry
!    nvtx,nvtxkm,ijkvtx,iwid,jwid,kwid,ncells,ijkcell ... mesh index descriptors
!     cdlt,sdlt ... cos,sin of toroidal angle subtended by the grid
!    c1x..c8z ... geoemtric coefficients for numerical differentiation
!     on a grid
!     on exit
!
!*** ROUTINES CALLED
!     gradc, divv, bcperv, axisgrad, mudote
!***END PROLOGUE AVEC
!
      cartesian = .T. 
!
!     calculate the curl of E at cell centers and store in curlc
!
!     AVEC calls CULRC, CURLV, and MUDOTE
!
      if (periodicx) then 
         call list (1, ibp1, 1, jbp1, 2, kbp1, iwid, jwid, kwid, nvtx, ijkvtx) 
      else 
         call list (2, ibp1, 1, jbp1, 2, kbp1, iwid, jwid, kwid, nvtx, ijkvtx) 
      endif 
!
!
!     EX COMPONENT
!
      call gradc (nvtx, ijkvtx, kbp1, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
         , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z&
         , c4z, c5z, c6z, c7z, c8z, vol, ex, exmu, eymu, ezmu) 
!
!     calculate the divergence of gradc at cell vertices and store in curlv
!
      call divv (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, exmu, eymu, ezmu, curlvx) 
!
!     ey component
!
      call gradc (nvtx, ijkvtx, kbp1, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
         , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z&
         , c4z, c5z, c6z, c7z, c8z, vol, ey, exmu, eymu, ezmu) 
!
!     calculate the divergence of gradc at cell vertices and store in curlv
!
      call divv (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, exmu, eymu, ezmu, curlvy) 
!
!
!     EZ COMPONENT
!
      call gradc (nvtx, ijkvtx, kbp1, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
         , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z&
         , c4z, c5z, c6z, c7z, c8z, vol, ez, exmu, eymu, ezmu) 
!
!     calculate the divergence of gradc at cell vertices and store in curlv
!
      call divv (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, exmu, eymu, ezmu, curlvz) 
!
      call list (2, ibp1 + 1, 2, jbp1 + 1, 2, kbp1 + 1, iwid, jwid, kwid, nvtx&
         , ijkvtx) 
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         curlvx, curlvy, curlvz, periodicx) 
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, curlvx, curlvy, curlvz) 
!
!     normalize axis divergence
!
      if (.not.cartesian) then 
         do j = 2, jbp1 + 1 
            do i = 2, ibp1 + 1 
!
               ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + kwid 
               rvvol = 1./vvol(ijk) 
               curlvx(ijk) = curlvx(ijk)*vvolaxis(ijk)*rvvol 
               curlvy(ijk) = curlvy(ijk)*vvolaxis(ijk)*rvvol 
               curlvz(ijk) = curlvz(ijk)*vvolaxis(ijk)*rvvol 
!
!
            end do 
         end do 
      endif 
!
!
!     project susceptibility tensor on to E
!
      call mudote (nvtx, ijkvtx, nsp, itdim, qdnv, qom, bxv, byv, bzv, dt, &
         transpos, clite, ex, ey, ez, exmu, eymu, ezmu) 
!
      dtsq = dt**2 
      dtc = clite*cntr*dt 
      dtcsq = dtc**2 
!test      DTCSQ=0.0
!
!
      if (cartesian) then 
!dir$ ivdep
         do n = 1, nvtx 
!
 
            ijk = ijkvtx(n) 
!
            hx(ijk) = (-dtcsq*curlvx(ijk)) + vvol(ijk)*(0.5*cntr*dtsq*exmu(ijk)&
               +ex(ijk)) 
            hy(ijk) = (-dtcsq*curlvy(ijk)) + vvol(ijk)*(0.5*cntr*dtsq*eymu(ijk)&
               +ey(ijk)) 
            hz(ijk) = (-dtcsq*curlvz(ijk)) + vvol(ijk)*(0.5*cntr*dtsq*ezmu(ijk)&
               +ez(ijk)) 
!
         end do 
!
      else 
!dir$ ivdep
         do n = 1, nvtx 
!
 
            ijk = ijkvtx(n) 
!
            hx(ijk) = (-dtcsq*curlvx(ijk)) + vvolaxis(ijk)*(0.5*cntr*dtsq*exmu(&
               ijk)+ex(ijk)) 
            hy(ijk) = (-dtcsq*curlvy(ijk)) + vvolaxis(ijk)*(0.5*cntr*dtsq*eymu(&
               ijk)+ey(ijk)) 
            hz(ijk) = (-dtcsq*curlvz(ijk)) + vvolaxis(ijk)*(0.5*cntr*dtsq*ezmu(&
               ijk)+ez(ijk)) 
!
         end do 
      endif 
!
!
      return  
      end subroutine avec 


      subroutine axisgrad(ibp1, jbp1, iwid, jwid, kwid, gradxf, gradyf, gradzf) 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ibp1 
      integer  :: jbp1 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      real  :: gradxf(*) 
      real  :: gradyf(*) 
      real  :: gradzf(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: j, i, ijk 
      real :: gradx, grady, gradz 
!-----------------------------------------------
!
!
!     calculate gradient at the axis by averaging gradients in the poloidal angle
!
!     WARNING!!! THIS VERSION ONLY WORKS FOR CARTESIAN GEOETRIES
      return  
!
      end subroutine axisgrad 


      subroutine axisvec(ibp1, jbp1, iwid, jwid, kwid, vvol, vvolaxis, vecx, &
         vecy, vecz) 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ibp1 
      integer  :: jbp1 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      real  :: vvol(*) 
      real  :: vvolaxis(*) 
      real  :: vecx(*) 
      real  :: vecy(*) 
      real  :: vecz(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: j, i, ijk 
      real :: avgx, avgy, avgz, rvvol 
!-----------------------------------------------
!
!
!     calculate average vector at the axis by volume averaging
!
!     WARNING!!! THIS VERSION ONLY WORKS FOR CARTESIAN GEOETRIES
      return  
!
      end subroutine axisvec 


      subroutine axisvol(ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, ijkvtx, vol&
         , vvolaxis) 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ibp1 
      integer , intent(in) :: jbp1 
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      integer , intent(in) :: nvtx 
      integer , intent(in) :: ijkvtx(*) 
      real , intent(in) :: vvol(*) 
      real , intent(in) :: vol(*) 
      real , intent(out) :: vvolaxis(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, j, i, ijk 
!-----------------------------------------------
!
!
      vvolaxis(ijkvtx(:nvtx)) = vvol(ijkvtx(:nvtx)) 
!
!
!     WARNING!!! THIS VERSION ONLY WORKS FOR CARTESIAN GEOETRIES
!
      do j = 2, jbp1 + 1 
!
         vvolaxis(iwid+1+(j-1)*jwid+kwid:ibp1*iwid+1+(j-1)*jwid+kwid:iwid) = &
            0.125*(vol(iwid+1+(j-1)*jwid+kwid:ibp1*iwid+1+(j-1)*jwid+kwid:iwid)&
            +vol(1+(j-1)*jwid+kwid:(ibp1-1)*iwid+1+(j-1)*jwid+kwid:iwid)+vol(1+&
            (j-2)*jwid+kwid:(ibp1-1)*iwid+1+(j-2)*jwid+kwid:iwid)+vol(1+iwid+(j&
            -2)*jwid+kwid:ibp1*iwid+1+(j-2)*jwid+kwid:iwid)) 
!
      end do 
!
      return  
      end subroutine axisvol 

